% Details for implementing normal-inverse-Wishart prior (see, for example,
% Arias, Rubio-Ramirez and Waggoner (2018)).

%% Place this section before the loop for the posterior sampler.

% Define prior parameters.
nuBar = 0; % Prior degrees of freedom
PhiBar = zeros(N); % Prior scale
PsiBar = zeros(M,N); % Prior location of regression coefficients
% Could add ones to PsiBar to represent a random-walk prior.
% Inverse of prior variance of regression coefficients (zero matrix implies
% diffuse prior).
OmegaBarInv = zeros(M); 

% Derive posterior parameters.
nuTilde = T + nuBar; % Degrees of freedom
OmegaTilde = (XX'*XX + OmegaBarInv)\eye(M); % Variance.
cholOmegaTilde = chol(OmegaTilde,'lower');
OmegaTildeInv  =  XX'*XX  + OmegaBarInv;
PsiTilde = OmegaTilde*(XX'*YY + OmegaBarInv*PsiBar); % Location
PhiTilde = YY'*YY + PhiBar + PsiBar'*OmegaBarInv*PsiBar - ...
    PsiTilde'*OmegaTildeInv*PsiTilde; % Scale
PhiTilde = (PhiTilde+PhiTilde')*0.5; % Ensure matrix symmetric

%% Replace the posterior sampler with this block.
% Draw reduced-form parameters from posterior.
phi.Sigma = iwishrnd(PhiTilde,nuTilde);
phi.Sigmatr = chol(phi.Sigma,'lower');
phi.Sigmatrinv = phi.Sigmatr\eye(N);
BBdraw = kron(phi.Sigmatr,cholOmegaTilde)*randn(M*N,1) + ...
    reshape(PsiTilde,N*M,1); 
phiDraw = phiDraw + 1;